POP Steering parameterization

Eddy diffusivity implementation by Bates, Tulloch, Marshall, and Ferrari

POP Steering Parameterization

Here are my notes on the development of a parameterization for eddy diffusivity. The scheme is based on the paper by Michael Bates, Ross Tulloch, John Marshall, AND Raffaele Ferrari. I'm working with Gokhan Danabasoglu and Matt Long here at NCAR and John Marshall at MIT.

Gokhan's Implementation Notes on Diffusivity with Steering Level Suppression

12 November 2014

Starting with $$K=u_{rms}∗L_{mix}$$

the general form for diffusivity, $K$, is given by equation (6) of Bates et al. (2014). Namely,

$$L_{mix} = {\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} $$

Here,

  • $u_{rms}$ is the root-mean-square (rms) eddy velocity;

  • $L_{mix}$ is the mixing length;

  • $L_{eddy}$ is the eddy diameter (depth independent);

  • $u_{mean}$ is the mean zonal velocity (resolved);

  • $c$ is the zonal eddy phase speed (depth independent);

  • $\Gamma = 0.35$;

  • $b1 \sim 4$.

Eddy Length Scales:

Rossby deformation radius = $L_r = {c_r \over |f|}$

Equatorial Rossby deformation radius $= L_{req} = {\sqrt{c_r \over 2\beta}}$

Rhines scale $= L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma_{vi} \over \beta}$

$c_r$ is the first baroclinic wave speed computed following equation (2.2) of Chelton et al. (1998) with $m=1$;

$f$ is the Coriolis parameter;

$\beta$ is the latitudinal variation of the Coriolis parameter; and

$\sigma_{vi}$ is the Eady growth rate given by

$\sigma_{vi} = {f \over \sqrt{R_i}}$

with $R_i$ the vertically integrated (over the 100 – 2000 m depth range) Richardson number.

So, any one these length scales could be used as an eddy length scale. An alternative is

$L_{eddy} = min (L_r, L_{req}, L_{Rh}).$

Eddy Velocity:

$u_{rms} = alpha*\sigma_{vi}*L_r$

where $\sigma_{vi}$ is the Eady growth rate based on local Richardson number and $\alpha$ is a scaling constant.

Zonal phase speed:

$c = - \beta * L_r^2$

References:

Bates, M., R. Tulloch, J. Marshall, and R. Ferrari, 2014: Rationalizing the spatial distribution of mesoscale eddy diffusivity in terms of mixing length theory. J. Phys. Oceanogr., 44, 1523-1540, doi: 10.1175/JPO-D-13-0130.1.

Chelton, D. B., R. A. deSzoeke, M. G. Schlax, K. E. Naggar, and N. Siwertz, 1998: Geographical variability of the first baroclinic Rossby radius of deformation. J. Phys. Oceanogr., 28, 433-460.

Tullock, R., J. Marshall, and K. S. Smith, 2009: Interpretation of the propagation of surface altimetric observations in terms of planetary waves and geostropic turbulence. J. Geophys. Res., 114, C02005, doi: 10.1029/2008JC005055.

Questions:

  1. In the 2D implementation, $u_{rms}$ and $u_{mean}$ specifications: upper-ocean vertically or integrated or at $z = 0$?

    A: $U_{rms}$ is not depth dependent; both $u_{rms}$ and $u_{mean}$ are for surface only

  2. In the 2D implementation, vertical profile will be specified by $N2(z)$?

    A: Yes, $N^2(z) \over N_{ref}(z)$ to be more precise

  3. Local $R_i$ use imbedded in sigma in $u_{rms}$ calculation?

    A: No, $u_{rms}$ is depth independent

  4. alpha = ?

    A: trial and error

  5. Cancellation of $f$’s in $u_{rms}$ calculation?

  6. Zonal phase speed equation correct? Both $\beta$ and $L_r$ will be positive, producing $c < 0$ always. This appears to be in contrast with Tullock et al. (2009).

    A: What is plotted in Bates is (U-c)

Parameterizing the Eddy Length Scale

$$K=u_{rms}∗{\Gamma * \color{red}{L_{eddy}} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} $$

The three length scales under consideration:

  1. Rossby Deformation Radius $\quad\quad L_r = {c_r \over |f|}$
  2. Rossby Equatorial Radius $ \quad\quad= L_{req} = {\sqrt{c_r \over 2\beta}}$
  3. Rhine's Scale $ \quad\quad L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma_{vi} \over \beta}$

The Rossby Deformation Radius $L_r$ depends on the the Baroclinic wave speed $c_r$ and the Coriolis force $f$. The CESM value of the first Baroclinic wave speed is derived as per eq 2.2 in Chelton 1998.

First Baroclinic Wave Speed $c_r$ Chelton vs CESM

In [1]:
%pylab inline
%run steering_setup ;
basedir='/scratch/jet'
###nc = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.007.pop.h.nday1.yr.0249.nc')
nc = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.007.pop.h.nday1.0249-01.nc')
omega=7.2921e-5
deg2rad=pi/180
f=2*omega*sin(lat*deg2rad)
Populating the interactive namespace from numpy and matplotlib

/loan/jet/install/anaconda/lib/python2.7/site-packages/mpl_toolkits/__init__.py:2: UserWarning: Module argparse was already imported from /loan/jet/install/anaconda/lib/python2.7/argparse.pyc, but /loan/jet/install/anaconda/lib/python2.7/site-packages is being added to sys.path
  __import__('pkg_resources').declare_namespace(__name__)

In [2]:
### get Variable from file c_rossby is in cm/s convert to m/s
c_rossby = nc.variables['C_ROSSBY'][0,:,:]


#Resample (aka re-project, re-grid) the NCEP data to target grid. First with nearest neighbour resampling...
c_rossby_nearest = pyresample.kd_tree.resample_nearest(orig_def, c_rossby, \
        targ_def, radius_of_influence=500000, fill_value=None)
c_rossby_mps = c_rossby_nearest/100.

fig1=MapZoneContour(lon_targcyc,lat_targcyc,c_rossby_mps,figsize=(14,9),
                    title="CESM First Baroclinic Wave Speed (m/s)",
                    fmt='%.1f' )
c_r_chelton=Image(filename=basedir+'/steering/chelton_sfig1.jpg')
display(c_r_chelton)
/loan/jet/install/anaconda/lib/python2.7/site-packages/matplotlib/text.py:52: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  if rotation in ('horizontal', None):
/loan/jet/install/anaconda/lib/python2.7/site-packages/matplotlib/text.py:54: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  elif rotation == 'vertical':

First Baroclinic Rossby Radius Chelton vs CESM

CESM Rossby deformation radius = $L_{eddy}$

$L_r = {c_r \over |f|}$

$L_{req} = {\sqrt{c_r \over 2\beta}}$

$\require{cancel} \cancel {L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma \over \beta}}$

$L_{eddy} = min (L_r, L_{req}, \cancel{L_{Rh}}).$

In [3]:
fcor = nc.variables['FCOR'][0,:,:]
fcor_nearest = pyresample.kd_tree.resample_nearest(orig_def, fcor, \
        targ_def, radius_of_influence=500000, fill_value=None)
absf=np.abs(fcor_nearest)
l_rossby_calc=c_rossby_nearest/absf

btp = nc.variables['BTP'][0,:,:]
btp_nearest = pyresample.kd_tree.resample_nearest(orig_def, btp, \
        targ_def, radius_of_influence=500000, fill_value=None)

l_rossbyeq_calc=np.sqrt(c_rossby_nearest/(2*btp_nearest))
l_eddy_calc=np.minimum(l_rossby_calc,l_rossbyeq_calc)

### convert regridded radius from cm to km
cm2km=1.e-5
l_eddy_calc_km=l_eddy_calc*cm2km
clevs=[10,20,30,40,50,60,80,100,150,230]
fig2=MapZoneContour(lon_targcyc,lat_targcyc,l_eddy_calc_km,figsize=(14,9),
                    levels=clevs,
                    title="CESM First Baroclinic Rossby Radius (km) $L_{eddy} \sim min(l_r,l_{req})$",
                    fmt='%.1i' )
l_r_chelton=Image(filename=basedir+'/steering/chelton_fig2.jpg')
display(l_r_chelton)
/loan/jet/install/anaconda/lib/python2.7/site-packages/numpy/ma/core.py:937: RuntimeWarning: overflow encountered in multiply
  result = self.f(da, db, *args, **kwargs)

Parameterizing the Zonal Eddy Phase Speed $\color{red}{c}$

$K=u_{rms}∗{\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - \color{red}{c}|^2 /u_{rms}^2 (z=0)} $

$\cancel{c = - \beta * {L_r^2}}$ $L_r$ too high at equator

$c = - \beta * L_{eddy}^2$

In [4]:
#beta
clevs=arange(1.e-14,25.e-14,1e-14)
norm=matplotlib.colors.Normalize(),
fig12=MapContourf(lon_targcyc,lat_targcyc,btp_nearest,
                    levels=clevs,
#                    norm=norm,
                    addzonal=True,
                    figsize=(14,5),
                    extend='both',
                    title="Beta (1/cm*s)")

clevs=np.logspace(1,15,30)
norm=LogNorm()
l_eddy_calc_sq=l_eddy_calc*l_eddy_calc
fig8=MapContourf(lon_targcyc,lat_targcyc,l_eddy_calc_sq,
                    addzonal=True,
                    levels=clevs,
                    norm=norm,
                    figsize=(14,5),
                    extend='both',
                    title="$L_{eddy}^2 (cm^2)$ log scale")

c=btp_nearest*l_eddy_calc*l_eddy_calc # units=1/(cm s) * cm^2 = cm/s

#clevs=np.linspace(-10,30,10)
clevs=arange(-10,30,1)
norm=LogNorm()
fig12=MapContourf(lon_targcyc,lat_targcyc,c,
                    addzonal=True,
                    levels=clevs,
#                    norm=norm,
                    figsize=(14,5),
                    contourlines=True,
                    extend='both',
                    title="eddy phase speed c = beta*leddy**2 (cm/s) linear scale(NOTE I'm plotting positive c here - changed to negative after this plot")

c=-1.*c
<string>:240: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.
/loan/jet/install/anaconda/lib/python2.7/site-packages/matplotlib/contour.py:1515: UserWarning: Log scale: values of z <= 0 have been masked
  warnings.warn('Log scale: values of z <= 0 have been masked')

Hughes Phase Speed (cm/s) from Tulloch Marshall Smith '09

In [5]:
hughes_c=Image(filename=basedir+'/steering/Hughes_phase_speed.png')
display(hughes_c)

Parameterizing the Zonal Velocity $\color{red}{u_{mean}}$

$K=u_{rms}∗{\Gamma * L_{eddy} \over (1 + b1 * |\color{red}{u_{mean}} - c|^2 /u_{rms}^2 (z=0)} $

In [6]:
#################  PLOT (u-c)^2  convert m^2/s^2 #######################
umean=nc.variables['U_MEAN'][0,:,:]
umean10_nearest = pyresample.kd_tree.resample_nearest(orig_def, umean,
                targ_def, radius_of_influence=500000, fill_value=None)

umean10_mps=umean10_nearest*1.e-2

clevs=np.linspace(-60,60,40)
fig11=MapContourf(lon_targcyc,lat_targcyc,umean10_nearest,
                    addzonal=True,
                    figsize=(14,5),
                    levels=clevs,
                    extend='both',
                    title="U zonal velocity (cm/s)")
pylab.xlim([-20.,20.])

# get 2D versions of the lat and lon variables add longitude start point here!
lon_ecco_orig, lat_ecco_orig = np.meshgrid(lon_ecco[:], lat_ecco[:])

orig_ecco_def = pyresample.geometry.GridDefinition(lons=lon_ecco_orig, lats=lat_ecco_orig)

ecco_u_nearest = pyresample.kd_tree.resample_nearest(orig_ecco_def, ecco_u, \
        targ_def, radius_of_influence=500000, fill_value=None)
ecco_u_nearest_masked=np.ma.masked_where(ecco_u_nearest<-10, ecco_u_nearest, copy=True)

#print ecco_u_cyc_masked[:,0]

clevs=np.linspace(-60,60,40)
#ecco_u_nearest = pyresample.kd_tree.resample_nearest(orig_def, ecco_u, \
#        targ_def, radius_of_influence=500000, fill_value=None)
fig11=MapContourf(lon_targcyc,lat_targcyc,ecco_u_nearest_masked*100,
                    addzonal=True,
                    figsize=(14,5),
                    levels=clevs,
#                    extend='both',
                    title="ECCO U zonal velocity (cm/s)")
pylab.xlim([-20.,20.])
Out[6]:
(-20.0, 20.0)
In [7]:
clevs=-np.logspace(3,-3,20)
print clevs
norm=matplotlib.colors.SymLogNorm(linthresh=1e-4, vmin=-1000, vmax=-.001)
fig12=MapContourf(lon_targcyc,lat_targcyc,c,
                    addzonal=True,
                    levels=clevs,
                    norm=norm,
                    figsize=(14,5),
#                    extend='both',
                    title="eddy phase speed c = beta*leddy**2 (cm/s) log scale")
pylab.xlim([-200.,0.])

uminusc=umean10_nearest-c  #units cm/s
#clevs=np.logspace(-2,3,10)
clevs=np.linspace(0,25,25)
#norm=LogNorm()
fig13=MapContourf(lon_targcyc,lat_targcyc,uminusc,
                    addzonal=True,
#                    norm=norm,
                    levels=clevs,
                    figsize=(14,5),
#                    extend='both',
                    title="U-c (cm/s) (max = %i/min = %i) (c calculated using Leddy)"%(uminusc.max(),uminusc.min()))
pylab.xlim([0,25])
[ -1.00000000e+03  -4.83293024e+02  -2.33572147e+02  -1.12883789e+02
  -5.45559478e+01  -2.63665090e+01  -1.27427499e+01  -6.15848211e+00
  -2.97635144e+00  -1.43844989e+00  -6.95192796e-01  -3.35981829e-01
  -1.62377674e-01  -7.84759970e-02  -3.79269019e-02  -1.83298071e-02
  -8.85866790e-03  -4.28133240e-03  -2.06913808e-03  -1.00000000e-03]

Out[7]:
(0, 25)
In [8]:
bates_uminusc=Image(filename=basedir+'/steering/bates2013_uminusc.png')
display(bates_uminusc)
In [9]:
uminusc_mps_sq=uminusc*uminusc*1.e-4   #units cm/s*cm/s*m/100cm*m/100cm
clevs=np.logspace(-4,1,40)
norm=LogNorm()
fig14=MapContourf(lon_targcyc,lat_targcyc,uminusc_mps_sq,
                    addzonal=True,
                    levels=clevs,
                    figsize=(14,5),
                    norm=norm,
                    #                    extend='both',
                    title=" ${(U-c)}^2 {(m^2/s^2)}$ max = %g/min = %g (log scale) "%(uminusc_mps_sq.max(),uminusc_mps_sq.min()))


clevs=np.linspace(0,.06,25)
fig14=MapContourf(lon_targcyc,lat_targcyc,uminusc_mps_sq,
                    addzonal=True,
                    levels=clevs,
                    figsize=(14,5),
#                    extend='both',
                    title=" ${(U-c)}^2 {(m^2/s^2)}$ max = %g/min = %g "%(uminusc_mps_sq.max(),uminusc_mps_sq.min()))
pylab.xlim([0,.06])
Out[9]:
(0, 0.06)
In [10]:
bates_uminuscsq=Image(filename=basedir+'/steering/Bates_u_minus_c_sq.jpg')
display(bates_uminuscsq)

Parameterizing the root-mean-square eddy velocity $\color{red}{u_{rms}}$ and the Eady Growth Rate $\color{red}{\sigma_{vi}}$

$K=\color{red}{u_{rms}}∗{\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} $

$u_{rms} = alpha*\color{red}{\sigma_{vi}}*L_{eddy}$

The original derivation of $\sigma_{vi}$:

$\sigma_{vi} = {f \over \sqrt{R_i}}$

with $R_i$ the vertically integrated (over the 100 – 2000 m depth range) Richardson number.

An alternate derivation of $\sigma_{vi}$ is:

$ R_i = {N^2 \over { ( \frac {\partial u}{\partial z} )^2 + ( \frac {\partial v}{\partial z} )^2} }$

$ N^2 = {{-g \over \rho_0 }\frac {\partial \rho}{\partial z}} $

After hydrostatic and geostrophic approximations

$f \frac {\partial v}{\partial z} = {{-g \over \rho_0 }\frac {\partial \rho}{\partial x}}; \quad\quad f \frac {\partial u}{\partial z} = {{g \over \rho_0 }\frac {\partial \rho}{\partial y}} $

so

$\frac {\partial v}{\partial z} = {{-1\over f}{g \over \rho_0 }\frac {\partial \rho}{\partial x}}; \quad\quad \frac {\partial u}{\partial z} = {{1\over f}{g \over \rho_0 }\frac {\partial \rho}{\partial y}} $

$\therefore$

$ R_i = {f^2 N^2 \over \underbrace{{ {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial y})^2 + {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial z})^2} }_\text{m^4}}\quad\quad = \quad\quad {f^2N^2 \over m^4} $

$\sigma_{vi} = {f \over \sqrt{R_i}}\quad = \quad {f \over \sqrt{{f^2N^2 \over m^4}}}\quad = \quad {{\cancel f m^2} \over \cancel f N}$

$RX_1 = RX_{east} = \Delta\rho_x = \rho_{i+1,j} - \rho_{i,j}$

$RY_1 = RY_{north} = \Delta\rho_y = \rho_{i,j+1} - \rho_{i,j}$

$RZ_1 = RZ_{k+1} = \Delta\rho_z = \rho_{k} - \rho_{k+1}$

$\displaystyle{1 \over L_{R_i}} \displaystyle\int_{2000m}^{100m} \left\lbrace { {-g\over\rho_0}{\frac {\partial \rho} {\partial z}} \over { {g^2 \over \rho_0^2 } \left[( \frac{\partial \rho}{\partial y})^2 + ( \frac{\partial \rho}{\partial z})^2\right] } } \right\rbrace dz$

Note: missing $f^2$ which will be cancelled when forming $\sigma_{vi}$

$\quad\quad$ so $\cdots$ this is not $R_i$

Implementation notes

For Level K

Numerator : Top $= -grav * RZ_{SAVE}(\cdots k+1) * dzwr(k)$

Denominator :

$ \begin{align} work1 = p25 & * ( RX(..,i_{east},k)^2 \\ & + RX(..,i_{west},k)^2 \\ & + RX(..,i_{east},k+1)^2 \\ & + RX(..,i_{west},k+1)^2 ) / DXT(i,j)^2 \\ \end{align} $

$ \begin{align} work2 = p25 & * ( RY(..,j_{north},k)^2 \\ & + RY(..,j_{south},k)^2 \\ & + RY(..,j_{north},k+1)^2 \\ & + RY(..,j_{south},k+1)^2 \\ \end{align} $

$work3 = {\left( TOP \over (grav^2*(work1+work2))\right)}*dzw(k)$

Notes:

1)Need to be careful at top and bottom of ocean
2)Accurate dzw(k) for each (i,j) to form $L_{R_i}$
3) When constructing $sigma$ itself, use $RZ_{SAVE}$ with a minimum N value
4) use eps2
In [11]:
sigma_old = nc.variables['SIGMA_AVG'][0,:,:]
sigma_old_nearest = pyresample.kd_tree.resample_nearest(orig_def, sigma_old, \
        targ_def, radius_of_influence=500000, fill_value=None)
clevs=arange(.005,.0625,.0025)

fig=MapContourf(lon_targcyc,lat_targcyc,sigma_old_nearest*86400,
                levels=clevs,
                figsize=(14,5),
                extend='both',
                title="Old Eady inverse time scale calc using alt sigma $(days^{-1}) \sim  {|f| \over \sqrt(R_i)}$")

sigma_avg1 = nc.variables['SIGMA_AVG1'][0,:,:]
sigma_avg1_nearest = pyresample.kd_tree.resample_nearest(orig_def, sigma_avg1, \
        targ_def, radius_of_influence=500000, fill_value=None)

clevs=arange(.005,.0625,.0025)
fig=MapContourf(lon_targcyc,lat_targcyc,sigma_avg1_nearest*86400,
                levels=clevs,
                figsize=(14,5),
                extend='both',
                title="New Eady inverse time scale calc using alt sigma $(days^{-1}) \sim  {{m^2} \over N}$")

tulloch_eady=Image(filename=basedir+'/steering/tulloch_eady.png')
display(tulloch_eady)
In [12]:
#################  PLOT u_rms with Leddy #######################
cm2m=.01
alpha=5.
u_rms_leddy = alpha*sigma_avg1_nearest*l_eddy_calc    # units = 1*1/s*cm = cm/s
u_rms_leddy_mps=u_rms_leddy*cm2m

clevs=np.linspace(0,.4,40)
fig99=MapContourf(lon_targcyc,lat_targcyc,u_rms_leddy_mps,
                    addzonal=True,
                    levels=clevs,
                    figsize=(14,5),
#                    extend='both',
                    title="$u_{rms} ({m \over s}) \quad alpha * \sigma * l_{eddy}$ (alpha=%i,max = %i/min = %i) "%(alpha,u_rms_leddy_mps.max(),u_rms_leddy_mps.min()))
pylab.xlim([0,.2])
Out[12]:
(0, 0.2)
In [13]:
bates_urms=Image(filename=basedir+'/steering/bates2013-urms.png')
display(bates_urms)
In [14]:
#################  PLOT u_rms^2  convert m^2/s^2 #######################
u_rms_leddy_mps_sq=u_rms_leddy_mps*u_rms_leddy_mps
#clevs=np.linspace(0,.06,50)
clevs=arange(0,.06,.001)
#clevs=np.logspace(-5,1,30)
norm=LogNorm()
fig99=MapContourf(lon_targcyc,lat_targcyc,u_rms_leddy_mps_sq,
                    addzonal=True,
                    levels=clevs,
#                    norm=norm,
                    figsize=(14,5),
#                    extend='both',
                    title="$u_{rms}^2 ({m^2 \over s^2})$ (alpha=%i,max = %g/min = %g) "%(alpha,u_rms_leddy_mps_sq.max(),u_rms_leddy_mps_sq.min()))
pylab.xlim([0,.06])
Out[14]:
(0, 0.06)
In [15]:
urms_bates=Image(filename=basedir+'/steering/Bates_urms_sq.jpg')
display(urms_bates)
In [16]:
uminusc_sq_ovr_urms_sq = uminusc_mps_sq/u_rms_leddy_mps_sq             #units cm/s /  cm/s  
uminusc_sq_ovr_urms_sq_masked=np.ma.masked_where(uminusc_sq_ovr_urms_sq >100., uminusc_sq_ovr_urms_sq, copy=True)
clevs=arange(0,5,.1)
fig14=MapContourf(lon_targcyc,lat_targcyc,uminusc_sq_ovr_urms_sq,                
                    addzonal=True,
                    levels=clevs,
#                    norm=norm,
                    figsize=(14,5),
                    extend='both',
#                    title="${(U-c)}^2 \over u_{rms}^2$ max = %g/min = %g "%(uminusc_sq_ovr_urms_sq.max(),uminusc_sq_ovr_urms_sq.min()))
                    title="${(U-c)}^2 \over u_{rms}^2$")
pylab.xlim([0,8])

clevs=np.logspace(-3,4,21)
norm=LogNorm()
fig14=MapContourf(lon_targcyc,lat_targcyc,uminusc_sq_ovr_urms_sq_masked,                
                    addzonal=True,
                    levels=clevs,
                    norm=norm,
                    figsize=(14,5),
                    extend='both',
                    title="${(U-c)}^2 \over u_{rms}^2$ (masking values over 100)")                    
#                    title="${(U-c)}^2 \over u_{rms}^2$ max = %g/min = %g "%(uminusc_sq_ovr_urms_sq.max(),uminusc_sq_ovr_urms_sq.min()))
pylab.xlim([0,100])

clevs=arange(0,5,.1)
fig14=MapContourf(lon_targcyc,lat_targcyc,uminusc_sq_ovr_urms_sq_masked,                
                    addzonal=True,
                    levels=clevs,
#                    norm=norm,
                    figsize=(14,5),
                    extend='both',
#                    title="${(U-c)}^2 \over u_{rms}^2$ max = %g/min = %g "%(uminusc_sq_ovr_urms_sq.max(),uminusc_sq_ovr_urms_sq.min()))
                    title="${(U-c)}^2 \over u_{rms}^2$ (masking values over 100)")
pylab.xlim([0,8])
Out[16]:
(0, 8)
In [17]:
bates_uminusc_sq_ovr_urms_sq =Image(filename=basedir+'/steering/Bates_suppress.jpg')
display(bates_uminusc_sq_ovr_urms_sq )

Parameterizing the Mixing term $\color{red}{L_{mix}}$

$\color{red}{L_{mix}} = \Gamma * L_{eddy} * Suppression$

$\Gamma = 0.35$

$Suppression= {1 \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2 )}$

$\color{red}{L_{mix}} = {\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)}$

In [18]:
gamma=.35
b1=4.0
supp=1./(1+b1*uminusc_sq_ovr_urms_sq)
clevs=np.linspace(0,1,30)
clevs=arange(0,1,.05)
norm=matplotlib.colors.Normalize()
fig=MapContourf(lon_targcyc,lat_targcyc,supp,
                    addzonal=True,
                    norm=norm,
                    levels=clevs,
                    figsize=(14,5),
                    title="Suppression factor ${1 \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2)}$")

pylab.xlim([0,.8])

l_eddy_calc_m=l_eddy_calc*.01

clevs=np.logspace(0,6,30)
norm=LogNorm()
fig8=MapContourf(lon_targcyc,lat_targcyc,l_eddy_calc_m,
                    addzonal=True,
                    levels=clevs,
                    norm=norm,
                    figsize=(14,5),
                    extend='both',
                    title="$L_{eddy} (m)$ log scale")


pylab.xlim([0,250e3])


lmix=gamma*l_eddy_calc_m*supp
clevs=np.logspace(0,5,40)
#clevs=np.linspace(0,2,30)#
#clevs=arange(0,20000,1000)
#norm=matplotlib.colors.Normalize()
norm=LogNorm()
fig=MapContourf(lon_targcyc,lat_targcyc,lmix,
                    addzonal=True,
                    norm=norm,
                    levels=clevs,
                    figsize=(14,5),
                    title="$L_{mix} = {{\Gamma * L_{eddy}} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2)}$")
#                    title="$L_{mix} = \Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2)}$")

pylab.xlim([0,10e3])
Out[18]:
(0, 10000.0)
In [19]:
bates_suppression =Image(filename=basedir+'/steering/bates2013_suppression.png')
display(bates_suppression )

Parameterizing Eddy diffusivity $\color{red}{K}$

$\color{red}{K}=u_{rms}∗L_{mix}$

In [20]:
kdiff=u_rms_leddy_mps*lmix
#clevs=np.linspace(0,1e3,30)

#clevs=arange(0,10000,500)
clevs=np.logspace(0,4,40)
norm=LogNorm()
fig=MapContourf(lon_targcyc,lat_targcyc,kdiff,
                    addzonal=True,
                    norm=norm,
                    levels=clevs,
                    figsize=(14,5),
                    title="Eddy diffusivity $K$")
In [21]:
bates_k =Image(filename=basedir+'/steering/bates2013_K.png')
display(bates_k )

Steering Level status

1. L_eddy values ok
2. Rossby wave speed ok
3. Eady inverse timescale ok
    a. Use alternate calculation otherwise has problems at equator and
    b.Richardson averages are too high causing sigma to drop too low
4. Zonal Eddy Phase Speed (c) 
    a. Too high 10S:10N
    b. Missing westward velocities in ACC which appear 
        in Hughes Data used in this parameterization
5. (U-c) Because of problem with c this term is also too high 10S:10N.
6. U_rms  (alpha*sigma*l_eddy)  OK        
    a. Using alpha of 5
7. Zonal mean velocity looks ok compared to ECCO annual average
8. Suppression OK
9. K is low by a factor of 5

Todo

1. Need to work on Zonal phase speed c
2. Plot vertical distribution of K
In [22]:
urms_bates2=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-002.jpg')
urms_bates3=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-003.jpg')
urms_bates4=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-004.jpg')
urms_bates5=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-005.jpg')
urms_bates6=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-006.jpg')
urms_bates7=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-007.jpg')
urms_bates8=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-008.jpg')
urms_bates9=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-009.jpg')
urms_bates10=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-010.jpg')
urms_bates11=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-011.jpg')
display(urms_bates3,urms_bates4,urms_bates5,urms_bates6,urms_bates7,urms_bates8,urms_bates9,urms_bates10,urms_bates11)

CESM $U_{mean} = U_{resolved}$ integraged over {10,50,100,200,500} meters

In [23]:
umean_cesm=Image(filename=basedir+'/steering/umean.jpg')
display(umean_cesm)
In [24]:
bates_len_vel=Image(filename=basedir+'/steering/bates_len_vel.jpg')
display(bates_len_vel)